Water at interface with proteins 



Giancarlo Pranzese 1 , Valentino Bianco 1 and Svilen Iskrov 1,2 

1 Departament de Fi'sica Fonamental, Universitat de Barcelona 
Diagonal 647, 08028 Barcelona, Spain 
2 Ecole Normale Superieure de Cachan, Paris 

E-mail: gfranzese@ub.edu 

Abstract. Water is essential for the activity of proteins. However, the effect of 
the properties of water on the behavior of proteins is only partially understood. 
Recently, several experiments have investigated the relation between the dynamics 
of the hydration water and the dynamics of protein. These works have generated a 
large amount of data whose interpretation is debated. New experiments measure the 
dynamics of water at low temperature on the surface of proteins, finding a qualitative 
change (crossover) that might be related to the slowing down and stop of the protein's 
activity (protein glass transition) , possibly relevant for the safe preservation of organic 
material at low temperature. To better understand the experimental data several 
scenarios have been discussed. Here, we review these experiments and discuss their 
interpretations in relation with the anomalous properties of water. We summarize the 
results for the thermodynamics and dynamics of supercooled water at an interface. We 
consider also the effect of water on protein stability, making a step in the direction 
of understanding, by means of Monte Carlo simulations and theoretical calculations, 
how the interplay of water cooperativity and hydrogen bonds interfacial strengthening 
affects the protein cold denaturation. 
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1. Introduction 

Water is ubiquitous is biological systems. It is a major component of cells and 
participates in the majority of the biological processes. It is usually considered essential 
for life, but it is still under debate why [T]. One possible reason is that water has many 
properties that are unusual with respect to other liquids [2J. 

The anomalous behavior of water is evident in the liquid phase. For example, 
fluctuations of volume and fluctuations of entropy have a minimum for liquid water, 
while in usual liquids they decrease when the temperature T is decreased. Volume 
fluctuations can be observed by measuring the compressibility K T , defined as how 
much the volume changes when the pressure P is changed at constant T, and entropy 
fluctuations are proportional to the specif heat Cp at constant P. For water at ambient 
pressure Kt has a minimum at 46° C and Cp at minimum at 35° C. 

The anomalies of water become more evident when T is decreased toward and 
below 0°C. For example, water has a maximum in density at 4°C Normal liquids, 
such as argon, reduce their density when the temperature decreases and reach their 
maximum density when they solidify in a crystal. Water, instead, below 4°C expands. 
Therefore, liquid water at 0°C has a density smaller than water at 4°C. By solidifying 
into ice, water expands even further, becoming "lighter". For this reason ice cubes float 
in a glass of water. This property has a dramatic consequence in processes such as 
the cryopreservation of biological cells, because the large amount of water in each cell 
expands when forms ice and breaks the cell. 

Another property of water that is relevant in these conditions is that water can 
stay in its liquid state even at T < 0°C, i.e. it can be supercooled below its melting 
temperature. Bulk water can be supercooled to — 41°C at atmospheric pressure, 
but in different conditions water can remain liquid even at lower temperatures. For 
example, it can be supercooled down to — 47°C when confined in vegetable fibers, or 
down to — 92°C when compressed at 2 kbars. An anti-freezing effect can be achieved 
also by dissolving in water polymers or proteins, with several practical applications. 
Understanding the details of this phenomena and how to regulate it could be very 
relevant for cryopreservation, food storage and refrigeration [31 H] . 

Also in its solid state, water is peculiar. Water is a polymorph with many different 
crystal phases, more than fifteen, some of which are stable only at pressure P > 100 GPa. 
But water can easily form a solid that is not a crystal. If quenched rapidly below 
— 123°C at ambient pressure, liquid water freezes in a metastable amorphous state, 
which is an arrested liquid configuration [2J. At low pressure water forms a low-density 
amorphous (LDA) state [5], while at high pressure it forms a high-density amorphous 
(HDA) state [6], separated by a volume discontinuity of ~ 27%, comparable to that 
between crystalline ice I and ice VI. A smaller discontinuity has been observed more 
recently [HE], but its interpretation is under debate. 

The discontinuities between amorphous states at different densities and the fact 
that quantities such as Kt or Cp largely increase in the supercooled state, show that 
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water has a complex behavior at low T. This observation and the fact that water remains 
in its liquid state at very low T when in contact with organic or inorganic interfaces 
suggest that water could play a main role in phenomena such as the so-called protein 
glass transition, or the protein cold denaturation at low temperature. 

1.1. Experiments 

To explore how the dynamics of proteins and water are related, S-H. Chen et al. in 2006 
studied by high-resolution quasi-elastic neutron scattering (QUENS) the structure and 
dynamics of water molecules in the hydration layer surrounding lysozyme proteins at 
temperatures around 220 K (— 53.15°C) [9]. Below this temperature the protein is in a 
solid-like "glassy state" , with no conformational flexibility and no biological functions. 
As the temperature is increased, the protein displays a harmonic atomic motion that, 
in hydrated proteins, suddenly becomes anharmonic and liquid-like at about 220 K. 
The change in the protein dynamics is believed to be triggered by the coupling with 
the hydration water through the hydrogen bonds because the hydration water displays 
a dynamic transition at a similar T [TUl fTT] . Chen interpreted the hydration water 
dynamic transition as a change in a main structural relaxation and recently extended 
this interpretation to hydrated proteins at high-pressure. 

This explanation has been questioned by Swenson et al. [12J. By using dielectric 
measurements on myoglobin in water-glycerol mixtures, they found a dynamic crossover 
at about 200 K and they interpreted it as an evidence that local (secondary) protein 
motions are controlled (slaved) by the local fluctuations in the hydration shell, as 
proposed by Fenimore et al. in 2004 [13J based on Mossbauer and neutron-scattering 
experiments. 

On the other hand, Pawlus et al. in 2008, based on measurement of conductivity 
on hydrated lysozyme, found no crossover around 220 K. They also ascribe the apparent 
crossover observed with QENS to a secondary relaxation and to a lack of resolution on 
main structural relaxation of QENS [H]. 

More recently nuclear magnetic resonance (NMR) experiments on water hydrating 
elastin and collagen, performed by Vogel, showed no crossover at 220 K but a crossover 
at about 200 K. This data have been interpreted as consistent with thermally activated 
tetrahedral jump motion of hydration water [15j. 

These and other experiments, therefore, show that it is difficult to achieve a clear 
understanding of the protein-hydration water coupling solely based on the experimental 
data. A possible way to gain further insight is offered by the numerical results of detailed 
computer simulations. 

1.2. Numerical results 

In 2006 Kumar et al. [16], by simulating lysozyme or DNA surrounded by water 
represented with the TIP5P model, found a dynamic transition of the macromolecules at 
about the same temperature of a dynamic crossover in the diffusivity of hydration water, 
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in the range 242 K-250 K. They also show numerical indications that this crossover 
coincides with the maximum of the isobaric specific heat of the whole system and the 
maximum fluctuations in tetrahedral order of the hydration water. 

This observation was confirmed in 2008 by Lagi et al. [T7J. They found a 
strong crossover in the water translational (a) relaxation time and in the inverse of 
its self-diffusion constant at about 223 K, by simulating hydrated lysozyme with water 
represented by the TIP4P-Ew model. They corroborate that the crossover corresponds 
to the maximum structural change, and they found a low activation energy, showing 
consistency with the neutron scattering and the NMR experiments. 

Nevertheless, in 2009 this view was disputed by Vogel [TS]. By simulating elastin or 
collagen hydrated by SPC-water he found only a weak crossover in the hydration water 
correlation time at about 200 K and a high activation energy associated to a secondary 
relaxation, consistent with activation energies determined in dielectric spectroscopy and 
NMR studies. 

It is evident from these and other numerical studies that the comparison of 
experiments with simulation is useful, but is not enough to elucidate the mechanisms 
that regulate the coupling of the water dynamics with the biomolecules dynamics. In 
particular, the difficulties in finding clear answers to the open questions rise for the fact 
that experiments and numerical results are both affected by errors and uncertainties. 
Therefore, it is natural to look for a theory that could be able to find exact relations 
and make predictions to test in further experiments. 

In the next section we describe a model for a hydration water monolayer that allows 
to develop a theory for these phenomena and to predict properties that could be verified 
in experiments. The model, moreover, allows for efficient simulations whose results 
complement the theoretical analysis. In section 3 we review results of this tractable 
model, including new data for the density and energy distributions at different pressures 
and very low temperature (subsection 3.3) and for the hydration percolation (subsection 
3.4). In section 4 we discuss new results about protein stability and confined water in 
the context of food processing and cells in living organism. In section 5 we give our 
conclusive remarks and we discuss possible extensions of this model. 

2. A tractable model for a hydration water monolayer 

We first consider the case of water nanoconfined between two hydrophobic surfaces. 
To fix the idea, let's consider the case with a distance 5 = 0.7 nm between the two 
surfaces. Because the confining surfaces are hydrophobic, the water molecules will not 
form hydrogen bonds (HB) with the surfaces. Experiments shows that in bulk each 
water molecule is surrounded by four nearest neighbor molecules at a distance of about 
3A with a structure that resembles that of a tetrahedron at low TV and P [19]. Hence, 
one would expect that each water molecule between the two hydrophobic plates will 
adjust in a way to form a distorted network of HBs with each molecules surrounded by 
other four. This has been indeed found by Kumar et al. [2D] by simulating TIP5P-water 
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in these conditions. 

To define a tractable model [2H [221 123J I2H I2S1 [2], we coarse-grain the structure 
described above, dividing the slab of space occupied by water into cells with volume 
v = 5r 2 , with a square section of size r > r , where r = 2.9A [26] is the closest 
approach (van der Waals) distance between water molecules (Fig. [I]). If the density of 
water is p, and the number of water molecules in N, the water volume is V = N/p. If 
the system is uniform and we divide the system in N cells, each cell has on average one 
water molecule. More in general if the system is not uniform, some cells can be empty 
or each cell can occupy a different volume Vi and the distance ry between two molecules 
in the cells j and j is the distance between their centers. For example, if i and j are 
the indices of nearest neighbor cells, the distance between the molecules in these cells is 
r l3 = (y^/S+y^J/S)/2 . 

We describe the isotropic (e.g., van der Waals) attractive and repulsive interactions 
between the molecules by a standard a Lennard- Jones interaction 



where e = 5.8 kJ/mol [27] is the attractive energy and the sum is over all the 
possible pairs of molecules % and j. Some modifications of this interaction, such as 
the introduction of a maximum cut-off distance or a hard-core distance, have been 
adopted in previous analysis of this model [2H EH ESI USB EH E21 [33] to simplify the 
numerical simulations. For the theoretical analysis, instead, this interaction term has 

been replaced by a more tractable discrete interaction [2Ill22ll2gll2HI25ll21l33]|351l33]. 

To take into account the energy and entropy variation when water molecules form 
HBs, we introduce for each water molecule i four bonding indices o~ij, one for each of 
the possible HBs with the four nearest neighbor water molecules j. The index can 
assume q different states, — 1, . . . , q, where q = 6 is a parameter whose value we will 
discuss in the following. Therefore, each water molecule has q 4 = 6 4 = 1296 possible 
states, that can be interpreted as possible rotational configurations. The total number 
of configurations for the system is q 4N , that for N = 10 5 (the maximum number that we 
have considered in our analysis) is an astronomical number (about 3 x 10 311260 ). When 
a molecule forms a HB the number of its accessible configurations decreases, and the 
energy of the system is reduced. 

To estimate the decrease in the number of accessible configurations, we observe that 
a HB is broken if it deviates from a linear bond more than ±30°. Therefore, only 1/6 of 
the whole continuous range of orientations [0, 360°] in the OH — O plane are associated to 
a bonded state. Hence, 5/6 of the possible configurations are not-bonded. By allowing 
q = 6 possible states for each bonding index cr^, we can count correctly the entropy loss 
associated to the formation of a HB if only one out of q states corresponds to a HB. 
This is achieved by allowing the formation of the HB between molecules i and j only if 
fi<rij,erji = 1j where by definition 8 a ^ = 1 if a = b and S a ^ = otherwise. 

Next, we take into account that a HB is broken if the OH — O distance is too large 
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[36. 37 1 . To simplify the analysis, we introduce a condition on the 0-0 distance r, 
allowing the formation of HB only if r < r max , with r max = Tqa/2 = 4.10A. Therefore, 
considering that the length of the OH covalent bond is toh = 0.96A [38J, the maximum 
bond length is r max — toh = 3.14A, consistent with other choices in literature [361 137] . 
In particular, we impose this condition by introducing a discrete variable for each cell 
i, with rii — 1 if rj < toa/2, otherwise rij — 0. A cell with rij = 1 is liquid-like because 
its (dimensionless) density pi = (Srl)/(5rf) > 1/2, while a cell with = is gas-like 
because has pi < 1/2. Hence, we define the number of HBs as 

N RB = J2 n i n A tJ ^ ( 2 ) 

(iJ) 

where the symbol denotes that the sum is performed over nearest neighbor cells 

i and j. 

From the experiments we know the formation of HBs leads to an open network 
of molecules with, on average, four neighbors instead of twelve as in argon-like fluids. 
The resulting volume per molecule with HBs is larger than the volume per molecule 
with no HBs. This is observable as the anomalous density decrease described in the 
introduction, that is a consequence of the formation of a macroscopic number of HBs. 
This effect is incorporated in the model by considering the total water volume to be 
given by 

V = Vo + iWHB, (3) 

where Vq is the water volume in absence of HBs, and ^hb is the increase of volume 
per HB. To estimate the parameter whb, we consider as reference values the increase 
between the density pih = 0.92 g/cm 3 of the ice Ih at atmospheric pressure and ice VI, 
with density pvi = 1-31 g/cm 3 , or ice ice VIII, with density pvm = 1-46 g/cm 3 . Ice Ih is 
characterized by hexagonal rings of HBs with an almost perfect tetrahedral structure, 
while ice VI and ice VIII have a structure consisting of two interpenetrating tetrahedral 
networks of HBs. The relative increase of density in these cases is equal to 0.42 and 
0.59, respectively. Hence, in the model we set the relative HB increase of volume per 
molecule equal to the average between the previous reference values, i.e. f hb / 5tq = 0.5. 

Note that this increase of volume per molecule is due to a decrease of first neighbors 
and does not imply an increase of distance between molecules. Our model, being coarse- 
grained, does not include all the details about the structure, but maintains the increase 
of volume per molecule with no effect on the distance r between molecules. In particular, 
the HB volume increase does not affect the calculation of Uo(r) of Eq.Q. 

The formation of a HB leads to an energy gain, represented in the model by an 
interaction term 

je HB = -JN UBl (4) 

where J is the characteristic energy of the covalent (directional) component of the HB. 
This term only accounts for the two-body component of the HB interaction. However, 
in water many-body effects are relevant and, in particular, the three-body term [39, 40J. 
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This can be observed from the T-dependence of the 0-0-0 angle distribution. This 
distribution becomes sharper around the tetrahedral angle when T decreases jH] . Hence, 
we include in the model the many-body (cooperative) effect due to HBs [121 B21 Hi], 
which minimizes the energy when the HBs of nearby molecules optimize the tetrahedral 
orientation. This is accomplished by further adding to the Hamiltonian in Eqs. ([!]) and 
Q the term 

=^coop = ~^ff ^ ^ n i } j 3<T ik ,a ie j (5) 
i {k,t)i 

where J a is the characteristic energy of the cooperative component of the H bond, and 
the sum over (k, £)j is performed over all the six different pairs of the four bonding 
indices of molecule i. 

As discussed by Stokely et al. in Ref . [33] , Eqs. Q and ^ have two free parameters: 
J and J a , respectively. Experiments estimate the HB in ice Ih to be ~ 3.0 kJ/mol 
stronger than in liquid water [15] . Attributing this increase to a cooperative interaction 
among HBs [IB], we can estimate the value of J a in the cell model to be ~ 1.0 kJ /mol, 
because for each HB there would be 6/2 pairs of J CT -interactions. The optimal HB 
energy, -Ehb ; has been measured to be ~ 23.3 kJ/mol [17] . By considering tetrahedral 
clusters of H-bonded molecules, with HB and van der Waals interactions up to the 
third nearest neighbor molecules, the value for the directional component of the HB 
is estimated by Stokely et al. |33J as J ~ 12.0 kJ/mol. Other experimental estimates 
suggest that breaking the directional component of the HB requires J ~ 6.3 kJ/mol [48J. 
It is, therefore, a reasonable estimate to set J a /J = 1/10 and to consider J as the only 
free parameter of the model. 

In the following we will briefly summarize some recent results for this model, to 
show that it reproduces in a qualitative way the properties of water. An appropriate 
choice of the free parameter J leads to results that, as for more detailed models, can be 
fairly rescaled on the known proprieties of a water monolayer. In this respect, this model 
is not better then detailed models, but, it has two features that detailed models have 
not. i) It is less computationally expensive, because it is coarse-grained. This allows 
to simulate very large number of water molecules (about a million) on simple desktop 
computers in a few hours, ii) More importantly and differently from the detailed models, 
this model is tractable for theoretical calculations. The trade-off for these advantages 
is that the model is coarse-grained and cannot give informations about some properties 
of water, such as the structure. Nevertheless, works is in progress to overcome this 
limitation. 

3. Results for water between hydrophobic plates 

We study the model described in the previous section by mean-field (MF) analysis and 
Monte Carlo (MC) simulations. The MF approach follows the Bethe-Peierls and the 
cavity method [49], by expressing the molar Gibbs free energy in terms of an exact 
partition function for a portion of the system, and taking into account the effect of all 



Water at interface with proteins 



8 



the rest of the system as a mean field acting on the border of this portion, as described 
in Ref. [21150]. 

MC simulations are performed at constant N, P, T, allowing the volume Vq in 
Eq.([3]) to fluctuate as a stochastic variable. To minimize the boundary effects, we 
consider periodic boundary conditions in the directions parallel to the confining surfaces. 
To study the thermodynamic properties of the model we adopt an efficient cluster MC 
dynamics, defined in Ref. [31] . or a continuous T algorithm, the histogram reweighting 
method, as in Ref. [21] . To study the dynamics of the HBs we adopt a standard 
Metropolis algorithm [35], while to study the diffusion properties we use the Kawasaki 
algorithm [2"9"] . 

3.1. The gas and liquid phases for nanoconfined water: transport properties 

For water nanoconfined between hydrophobic plates, the model display a gas-liquid 
first-order phase transition ending in a critical point C, that qualitatively resembles the 
gas-liquid transition for bulk water (Fig. [2^,). In their work, de los Santos et al. [301 [51] 
verify that the water diffusion constant D decreases when water changes from gas to 
liquid. The change in D is strong far from the critical point C, and it disappears at 
C. In the liquid phase, de los Santos et al. j3U] observe that, as for bulk water, the 
nanoconfined system has a region in the P-T plane where D increases for increasing P. 
This behavior is an anomaly of water, because in normal fluids D decreases for increasing 
P [521 [53] . This anomalous behavior is qualitatively rationalized as a consequence of 
the HB formation and it has been observed both in bulk and confined water [51]. 

At lower temperature, by decreasing T at constant P the model displays the line 
of temperatures of maximum density (TMD) j35j EH], as in the bulk case (Fig. [2^,). The 
TMD line in the P-T plane has a positive slope at low P and negative slope at high P. 

The nanoconfinement does not allows the formation of crystal ice. Nevertheless, it 
is possible to calculate where in the P-T plane D(P,T) is constant and to show that 
the lines of constant-/} qualitatively resemble the melting line of bulk water [301 EI] ■ 

3.2. Water monolayer compared to protein hydration water 

At low T the diffusion constant largely decreases and, eventually, the monolayer becomes 
subdiffusive, i.e. the mean square displacement of water molecules never reach the 
diffusive regime [301 EE] • This property has been observed experimentally below 320 K 
by neutron scattering in a monolayer of water hydrating a myoglobin surface at low 
hydration level (h = 0.35 g H^O/g of protein), corresponding to a number of water 
molecules sufficient to cover the entire protein surface [53] . 

Franzese and de los Santos in 2009 [29] found that a water monolayer partially 
hydrating a hydrophobic surface, described by the model considered here, would display 
a very slow dynamics for the HBs at low P and T. In these conditions the HB correlation 
function C(t), that quantifies how much the HBs are correlated in time, is almost not 
changing in time, showing that the water dynamics is completely frozen. Hence, water 
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is in its glassy state at low T and low P, consistent with the observed freezing of the 
incoherent intermediate scattering function of water hydrating myoglobin at hydration 
level h = 0.34 g H^O/g of protein at about 180 K [56]. Glassy water is observed also in 
bulk, but below 150 K. 

This slowing down of the dynamics is well understood in our model where, at low T 
and low P, the number of HBs largely increases when the T decreases. This progressive 
building up of the HB network traps the water molecules in a percolating network of 
HBs that leaves small areas of the hydrophobic surface completely dehydrated. 

At higher P the decreased number of HB allows to the small dry cavities to 
slowly equilibrate, due to large-scale rearrangements of the HBs and leading to partial 
dehydration of the surface. In this case water can slowly flow on the surface because 
the high pressure reduces the volume per molecules and disfavors the formation of 
HBs. Hence, the HB network builds up at lower T with respect to the low pressure 
condition [29]. In this case the long time behavior of C(t) is well described by a stretched 
exponential function 



where Co, r and /3 < 1 are fitting constant. For (3 = 1 the function is exponential, and 
the more stretched is the function, the smaller is the exponent < (3 < 1. 

Franzese and de los Santos [29J predicted that at low P, by increasing P, (i) the 
time needed for the HBs to decorrelate decreases, i.e. water can relax more rapidly, and 
(ii) the (3 exponent decreases going from (3 = 0.8 to (3 = 0.4 for a pressure approaching 
a characteristic value Pc- We will discuss further about Pc in the following. By 
increasing the pressure even further, for P > Pc>, the HB correlation function relaxes 
faster and the exponent becomes = 1, i.e. C(t) becomes an exponential function. 

The prediction about (3 is consistent with the experimental findings of Settles 
and Doster showing that, for water hydrating myoglobin at low hydration level, the 
incoherent intermediate scattering function at large Q vector, i.e. the dynamics of 
density fluctuations at short distance, are well described by stretched exponential 
functions with (3 varying between 0.4 and 0.3 at 320 K [55]. 

We observe that the theoretical lower limit for (3 is expected to be 1/3 [57J [58] and 
that 1 — (3 is a measure of the heterogeneity in the system. Therefore, the prediction [25] 
that the stretching parameter [3 approaches its smallest possible value when the pressure 
tends to Pc implies that the water monolayer reaches its maximum in heterogeneity at 
Pc- In the following we will discuss further this point. 

3.3. Thermodynamics at very low temperature 

To understand the origin of the heterogeneity at Pc, we recall here that the properties 
of water are consistent with theories that propose different mechanisms and different 
phase behaviors at very low temperature (approximately 150 K< T < 200 K). At these 
temperatures, supercooled bulk water forms ice, while confined water in appropriate 



C(t) = C exp -{t/r) 



(6) 
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conditions can be kept liquid [59). The different theories can be summarized in four 
possible scenarios for the P — T phase diagram. 

(i) In the stability limit (SL) scenario [60] it is hypothesized that the limits of 
stability of superheated-stretched liquid water changes its slope in the P-T plane from 
positive at high T, to negative at low T and negative P, giving as a consequence the 
anomalous increase at low T of quantities such as Kt, Cp and ap. 

(ii) In the liquid-liquid critical point (LLCP) scenario [HI] it is hypothesized the 
existence of a first-order phase transition line in the supercooled liquid region, with 
negative slope in the P-T plane and terminating in a critical point C. This phase 
transition separates two liquid phases, both metastable with respect to the crystal 
phases: one with low density, resembling the LDA disordered ice, and one with high 
density, resembling the HDA disordered ice. The low-density-liquid (LDL), high-density- 
liquid (HDL) critical point C has been predicted at positive P )61) or negative P [62), 
depending on the water model adopted in simulations. The anomalies of water are the 
consequence of approaching C . 

(iii) In the singularity-free (SF) scenario [63] it is hypothesized that the HBs have 
no cooperativity. In this case it is shown that the anomalous increase of Kt-, Cp and 
«p is a consequence of the low-T anticorrelation between volume and entropy, also 
responsible for negative slope in the P-T plane of the line of temperatures of maximum 
density (TMD). 

(iv) In the critical-point free (CPF) scenario [61] it is hypothesized that the LDL- 
HDL first-order phase transition line extends to P < 0, reaching the superheated limit 
of stability of liquid water and with no critical point. As a consequence the HDL has a 
superheated-stretched limit of stability similar to that predicted in scenario (i). 

As discussed by Stokely et al. [33J, it is still unclear which of the scenarios best 
describes water, because there is no definitive experimental test. In Ref. [33], the same 
tractable model presented in the previous section is analyzed by means of theoretical 
calculations and numerical simulations, showing that the four scenarios (i)-(iv) may 
be mapped in the space of the parameters J and J a , representing the strength of 
the HB directional component and the strength of the HB cooperative component, 
respectively. The relation J a /J = 1/10 discussed at the end of the previous section, 
and based on estimates from experimental data, supports the prediction of a liquid- 
liquid critical point C at positive pressure for supercooled water (Figj^). The model 
allows to distinguish the LDL phase and the HDL phase by their different densities 
and characteristic energies. The presence of two maxima in the distributions of these 
quantities marks the coexistence of the two phases and corresponds to the occurrence 
of two minima in the free energy of the system (Figj3|. 

At the critical point C the cooperativity of the HBs is maximum. Hence, 
cooperative rearrangements of the HB are necessary to allow the relaxation of the 
dynamics. These rearrangements occur on different length-scales, each associated with 
a different time-scale. As a consequence a single time-scale cannot be defined for the 
dynamics, resulting in a stretched decay of the correlation function, i.e. in heterogeneous 
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dynamics. It is, therefore, C the origin of the heterogeneity described in the previous 
section and occurring at Pc> of the liquid-liquid critical point. 

3.4- Hydration percolation 

It possible to study the cooperative regions, and their length-scales, by using a 
geometrical approach based on the concept of correlated percolation [HSl ,31J - We define 
a cluster of correlated water molecules as described in the following steps. 

• The first step consists in including in the cluster one of the bonding indices cry of 
a randomly selected water molecule i in the hydration shell. 

• The second step is to add to the cluster another bonding index of the same molecule 
i with probability p same = min {1, 1 — exp[— J a j (fegT)]} (where ks is the Boltzmann 
constant), or the facing bonding index cx^ of the nearest neighbor molecule j with 
probability pacing = min{l, 1 — exp[— J'/(/cbT)]} where J' = J — Pvhb- In this 
expression, J' is the P-dependent effective coupling between two facing indices as 
results from the enthalpy U + J^hb + ^coo P + PV of the system Eq.s ([T|-([5]). The 
quantity J' can be positive or negative depending on P. If J' > 0, an index can be 
added, with probability Pf ac ing, to the cluster only if it is in the same state as the 
other indices already in the cluster. Instead, if J' < 0, the index can be added only 
if it is in a different state with respect to the index to whom it will be connected. 

• The third step is to randomly select a bonding indices on the border of the cluster 
and pick at random one of the indices on the same molecule, or the facing index 
on a bonded molecule, that is not already in the cluster, and to include it in the 
cluster with probability p samc or Pf ac i ng , respectively. This step is repeated until all 
the possible directions of growth for the cluster have been considered. 

The resulting cluster statistically represents the region of correlated HBs, as can be 
shown [52], and its characteristic linear size statistically corresponds to the correlation 
length of the water molecules. Therefore, in the vicinity of the liquid-liquid critical point, 
where the correlation length increases, it is possible to observe that the size of the clusters 
of correlated water molecules increases. At the liquid-liquid critical point the correlation 
length diverges and a cluster of correlated water molecules spans (percolate) the whole 
monolayer. General results on correlated percolation theory allow to find mathematical 
relations between thermodynamics quantities and percolation quantities. In particular, 
it can been shown that the mean size of the clusters defined above diverges with the same 
power-law exponent as the compressibility of water and that the distribution of number 
n(s) of finite cluster of size s per water molecule is exponential far from the liquid-liquid 
critical point C, while follows a power law with exponent r near C From general 
considerations it is possible to show that r = 1 + d/Dp where d = 2 is the effective 
dimensionality of the monolayer and Dp is the fractal dimension of the clusters |67j. 
Preliminary estimate of r ~ 2 suggests that the clusters of correlated water molecules 
are compact with Dp ~ 2 (Fig Eh [UH]- Since the compressibility is proportional to 
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the density fluctuations, the clusters allows for a geometrical analysis of the diverging 
density fluctuations near the liquid-liquid critical point [32] . 

3.5. Hydrogen bonds dynamics on Hydrated Protein 

The density fluctuations are observable also far from the liquid-liquid critical point C . 
In particular, they can be observed along a line in the P-T phase diagram that emanates 
from C into the one-phase region and marks the maxima of the correlation length. This 
line has been named after Widom [691 121 EQ] and can be characterized in the study of 
hydrated protein. For example, Kumar et al. [28] used the tractable model described 
above to investigate the case of a percolating monolayer of water molecules adsorbed on 
a protein surface, with hydration level about h ~ 0.4 g H 2 0/g dry protein. Under these 
conditions the protein is immobile and inhibits the ice crystallization because it forces 
the water molecules out of the positions corresponding to a crystal configurations. The 
authors studied the HB dynamics, regardless if the HBs are formed within the water 
molecules or with the surface. 

They first locate the Widom line, by Monte Carlo simulations, and observed that it 
corresponds to the locus where there is the largest change in the number of HBs. At P 
and T above the Widom line, water has a few HBs, while at P and T below the Widom 
line, it has a well-developed network of HBs. This change of structure is reflected by 
a maxima in constant-P the specific heat along the Widom line. The simulations also 
reveal the presence of a dynamic crossover in the HB correlation function C (t) when the 
Widom line is crossed at constant P. Using mean field theory and making the hypothesis 
that the dynamics is dominated by the rearrangements of the HBs, the authors calculate 
the activation energy for the relaxation of the system and show that it gives the same 
relaxation time calculated by Monte Carlo simulations [35]. The proposed mechanism 
consists in breaking a HB that does not fit into the tetrahedral arrangement and reorient 
the molecule to optimize locally the tetrahedral configuration. This mechanism has been 
confirmed also by simulations of the hydration shells of elastin-like and collagen-like 
peptides [18]. In particular, Kumar et al. |35j predict (i) how this barrier is affected by 
the variation of P, (ii) how the crossover T is affected by the variation of P, and (iii) 
that for any P the correlation time at the crossover T is the same (isochronic crossover). 
Experiments with hydrated lysozyme, spanning a range of pressures going from ambient 
pressure up to 1600 bar, performed by the group of S.-H. Chen at MIT, have confirmed 
these three predictions [TTJ [31] . 

More recent analysis of lysozyme proteins at a lower hydration level (h =0.3 g 
H20/g dry protein) reveals another surprising results [72]. At this very low hydration, 
dielectric spectroscopy, probing the proton relaxation, displays that at ambient P not 
only there is at about 250 K the dynamic crossover described above, but also another 
crossover at about 180 K. The study of the tractable model presented here associates 
this lower-T crossover to the saturation of the cooperative ordering of the HB network. 
Specifically, the HBs rearrange to maximize the number of tetrahedral orientations 
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among the bonds. This ordering is marked by a new specific heat maxima that can be 
calculated in the model at about 180 K. Therefore, summarizing, the model predicts a 
broad specific heat maximum at about 250 K, due to the saturation of a macroscopic HB 
network, and a sharper specific heat maximum at about 180 K, due to the saturation of 
the tetrahedral ordering of the HBs. By increasing the pressure of the hydrated protein, 
Mazza et al. predict that these two maxima merge and then diverge at the liquid-liquid 
critical point C [75] . 

4. Discussion: Implications in food science 

The presence of a liquid-liquid critical point C at low T in confined water could be 
an undesirable property for the storage of frozen food and, more in general, biological 
cells. This is because in the vicinity of a critical point between two liquid phases, both 
metastable with respect to the crystal phase, large density fluctuations occur. These 
enhanced fluctuations would drastically changes the pathway for the formation of a 
crystal nucleus, because the crystal would form from the dense fluid, instead that from 
the low density fluid. As a consequence, there would be a strong reduction of the crystal 
nucleation free-energy barrier and, hence, an increase by many orders of magnitude of 
the crystallization rate, as theoretically predicted by ten Wolde and Frenkel [74] . 

Under this conditions the enhanced formation of ice could destroy the cells as a 
consequence of the increase of volume of ice with respect to the liquid water j3j 0]. 
Hence, the best way to preserve the cells would be to freeze them at a T that is far away 
from the liquid-liquid critical temperature. Nevertheless, further analysis show that 
the situation is even more complex. Indeed, simulations of a model with a metastable 
liquid-liquid critical point display enhanced crystallization rate not only in the vicinity 
of the critical point, but in the vicinity of the whole region of liquid-liquid coexistence 
[75J and possibly also in the one-phase region above the critical point along the Widom 
line (2]. Work is in progress to elucidate these implications. 

In the evaluation of these effects is extremely important to properly include the 
interaction with the confining surfaces. We are presently studying how to incorporate 
these effects in our tractable model. A first step in this direction is to include a 
description of the hydrophobic effect. Frank and Evans [76 J and Silverstein et al. 
[77] proposed that supercooled water forms highly structured "ice-like" regions in the 
hydration shell of nonpolar solutes. Stillinger [78] proposed that HBs in the hydration 
shell are not significantly perturbed near small hydrophobic solutes, while the HB 
network is strongly affected by hydrophobic particles with size above a characteristic 
value. Chandler estimated this value of the order of 1 nm on the basis of free 
energy calculations [79]. Muller explained the vibrational and NMR spectroscopy 
results by suggesting enthalpic strengthening of the hydration HBs with a simultaneous 
entropy increase in the hydration shell [HQ]. We are presently including the enthalpic 
strengthening in our model, and properly accounting for the entropy increase for the 
study of water in confined by hydrophobic nanoparticles. Our preliminary results [HI] 
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show a surprising change of thermodynamic fluctuations at low T, whose implication is 
a large decrease of compressibility also at very low nanoparticles concentration. This 
finding suggests that adding hydrophobic particles at low concentration in organic 
solutions would decrease the density fluctuations and the formation of ice. 

Further confirmation of the validity of our assumption about the hydrophobic effect 
comes from another study that we are performing to establish if our tractable model 
is able to describe the stability of proteins with respect to changes of temperature 
and pressure. In this study we consider how the hydrophobic interaction with water 
of a coarse-grained protein induces hot denaturation, folding, cold denaturation and 
pressure denaturation. Our preliminary results |82j display a region of stable folded 
configurations that, in the P-T phase diagram, has the same qualitative features of the 
experimental stability diagram of myoglobin [83J (Fig{5]). 

5. Conclusions 

We have introduced a tractable model for a water monolayer hydrating surfaces of 
proteins and, more generally, confined water. The model includes HB cooperativity 
and elucidates how the many-body component of the HB is important to understand 
the low-T behavior of water. In particular, parameters estimated from the experiments 
suggest the occurrence of a critical point at T ~ 180 K and P ~ 0.13 GPa at the end 
of a first order coexistence line between two liquids with different densities [231 EH [33] . 

The model shows that the liquid-liquid critical point affects the low-T dehydration 
of a hydrophobic surface [2S]. The cooperativity of water induces dynamic 
heterogeneities that reach their maximum when the cluster of correlated HBs percolate 
EH El EH]- 

This heterogeneous dynamic behavior is revealed by a strongly non-exponential 
relaxation of the HB dynamics and by a subdiffusive translational motion of the water 
molecules in the hydration shell [29], as observed in hydrated proteins at low T [55, 56J. 

The model predicts that, at low protein hydration, shell water should be 
characterized by two structural transitions. One associated to the macroscopic formation 
of HBs [351 Ell EH], occurring at about 250 K for low-hydrated lysozyme, and another 
associated to the tetrahedral reordering of the HBs, at about 180 K [72]. These two 
structural changes are at the origin of two dynamic crossovers: the one at higher T has 
been observed by QENS experiments on hydrated lysozyme [TTJ [31] , and both have been 
measured by dielectric spectroscopy on hydrated lysozyme [72] • The pressure behavior 
of these crossover is consistent with the presence of the liquid-liquid critical point C at 
high P and low T [HI [75]. 

The large increase of density fluctuations in the vicinity of the critical point C is 
expected to enhance the water crystallization process. As a consequence of the expansion 
of ice, this process could destroy biological structures in a crowed environment, as 
for example in food stored at low temperature [3j H]. It is, therefore, relevant to 
understand if this process is affected by confinement or if it could be controlled. Our 
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preliminary results [Hi] show that hydrophobic confinement has a strong effect of the 
thermodynamics of water, suppressing the density fluctuations associated to the critical 
point C This result suggests that by dissolving hydrophobic nanoparticles at low 
concentration it could be possible to control the water compressibility and the formation 
of large crystals. 

These predictions are based on a modelization of the hydrophobic interaction that 
is able to reproduce a stability diagram for a coarse-grained protein. In particular, the 
model shows that the protein cold denaturation and the pressure denaturation can be 
explained as a consequence of the strengthening of interfacial water-water HBs at the 
hydrophobic interface and the proper account of the entropy change due to the presence 
of the interface [82J. 

Work is in progress to include other features in the model and to use it to make other 
predictions in different contexts. For example, we are extending the model in such a 
way to describe the ice formation and analyze how the crystallization process is affected 
by interfaces. We are also developing the generalization in three dimensions to allow the 
study of many layers of hydration water and to extend the investigation to the bulk case. 
In this way our research about protein stability will be developed also in bulk water. 
All together these generalizations will allow us to explore situations of possible interest 
in food science and biology, as for example the effect on water structure and dynamics 
of preservative agents such as trehalose, or antifreeze proteins, or cryoprotectants such 
as glycerol or dimethyl sulphoxide. 
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Figure 1. Schematic representation of a water monolayer. Top view of water molecules 
with an oxygen atom (red) and two hydrogen atoms (blue) distributed over a surface, 
between two hydrophobic plates (not represented). Possible hydrogen bonds are 
represented by gray sticks. The total surface area is divided in equal-size square cells 
(dashed lines). In the tractable model adopted here, the coordinates of each molecule 
inside a cell are coarse-grained. A configuration of water molecules is represented by 
the occupancy state of the cell (local density) and the states of the four bonding indices 
of each molecule, accounting for the hydrogen bonds formed with water molecules in 
the nearest neighbor cells. 
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Figure 2. Density at constant pressure for a water monolayer nanoconfined between 
hydrophobic slabs of infinite section and separated by a distance S = 0.7 nm. Results 
are from Monte Carlo simulations of a system with N = 15625 water molecules. The 
pressure P is calculated in GPa, the temperature T is rescaled with respect to the 
liquid-gas critical point temperature Tc of the water model, the density p — N/V is 
rescaled by the volume vq = r^S ~ 0.059 nm 3 , where tq — 2.9A is the van der Waals 
distance, (a) Isobars at high T, around the liquid-gas critical point C, marked by a 
full large circle. At P < Pc — 0.18 GPa a discontinuity in the isobars denotes the 
coexistence of the gas (at pvo < 0.5) with the liquid (at pvo > 0.5). In the circled 
region for 0.02 GPa< P < 0.06 GPa the isobars reach the temperature of maximum 
density (TMD) at pv Q ~ 0.65 and T/T c ~ 0.35. (b) At much lower T, the isobars 
display another discontinuity in density that is evident at P = 0.18 GPa, and smoothly 
disappears, within the calculation error, for P < 0.13 GPa. This discontinuity denotes 
the coexistence of a liquid at higher density with a liquid at lower density. Where 
the discontinuity disappears, the system displays a liquid-liquid critical point C (large 
open circle, with Pc — 0.13 GPa and Tc> — 180 K. In both panels errors are of the 
order of the size of symbols. 




Figure 3. The histogram of rescaled density pvo and rcscalcd energy [U + ■JZhb + 
^coop]A above [panels (a)-(b)] and below [panels (c)-(d)] the liquid-liquid critical point 
C . In the one-phase region, density histogram (a) and energy histogram (b) display 
one single maximum. At the liquid-liquid phase separation, density histogram (c) and 
energy histogram (d) display two maxima separated by a minimum, corresponding 
to two coexisting phases with different densities and energies. Calculations are 
from Monte Carlo simulations of a monolayer with N — 15625 water molecules at 
T = 175.1 K and P = 0.12 GPa (a)-(b), and T = 173.8 K and P = 0.13 GPa (c)-(d). 
The liquid- liquid critical point C is estimated at Tc> — 180 K and Pc> — 0.13 GPa. 
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Figure 4. Histograms of number n(s) of clusters of size s near the liquid-liquid 
critical point C for a monolayer with N = 25600 water molecules. For T < Tc> (a) 
and T > Tc (c) the distribution of n(s) is exponential, while for T ~ Tc (b) it 
becomes approximately a power law, whose leading term is n(s) ~ s T with an estimate 
r ~ 2. Lines are fits of the histogram with exponential functions in (a) and (c) and 
with n(s) ~ s T in (b). 
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Figure 5. Folding and unfolding of a coarse-grained protein suspended in water at 
different temperatures T and high pressures P. Typical configurations in a system in 
two dimensions are represented for the sake of schematic description of the process. 
The protein is represented as a fully hydrophobic chain (in white), surrounded by 
water molecules (turquoise background), (a) At high pressure and high T, the protein 
is unfolded and the surrounding water has only a few hydrogen bonds (HBs) formed, 
represented as colored sticks. Different colors of the HBs correspond to different relative 
orientations of the HBs. (b) At the same pressure but lower T, the protein start to 
fold in a molten globule state, (c) At lower T the protein folds, while the surrounded 
water has a large number of HBs. (d) At much lower T we observe cold denaturation 
of the protein when the number of water HBs is largely reduced due to the combined 
action of T and P. (e) At higher P the denaturation effect is observed at higher T, 
and we observe the protein only in open configurations. 
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